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Abstract 

The nonlinear climbing sine map is a nonhyperbolic dynamical system exhibiting 
both normal and anomalous diffusion under variation of a control parameter. We 
show that on a suitable coarse scale this map generates an oscillating parameter- 
dependent diffusion coefficient, similarly to hyperbolic maps, whose asymptotic 
functional form can be understood in terms of simple random walk approximations. 
On finer scales we find fractal hierarchies of normal and anomalous diffusive regions 
as functions of the control parameter. By using a Green-Kubo formula for diffusion 
the origin of these different regions is systematically traced back to strong dynami- 
cal correlations. Starting from the equations of motion of the map these correlations 
are formulated in terms of fractal generalized Takagi functions obeying generalized 
de Rham-type functional recursion relations. We finally analyze the measure of the 
normal and anomalous diffusive regions in the parameter space showing that in 
both cases it is positive, and that for normal diffusion it increases by increasing the 
parameter value. 

Key words: Deterministic diffusion, nonhyperbolic maps, fractal diffusion 
coefficient, anomalous diffusion, periodic windows 
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1 Introduction 



Low- dimensional time-discrete maps are among the most important models for 
exploring different aspects of chaos. These systems display a very rich dynam- 
ical behavior but are still very amenable to straightforward computer simula- 
tions. Even more, in some cases rigorous analytical solutions are possible. After 
it was realized that diffusion processes can be generated by microscopic deter- 
ministic chaos in the equations of motion, time-discrete maps became useful 
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tools in deterministic transport theory. The analysis of these simple models re- 
quired to suitably combine nonequilibrium statistical mechanics with dynam- 
ical systems theory leading to a more profound understanding of transport in 
nonequilibrium situations [1,2,3,4]. However, time-discrete maps provide not 
only a suitable starting point for studying normal diffusion but also for investi- 
gating the anomalous case [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. 
Moreover, there are certain classes of more realistic models which share spe- 
cific properties of maps such as being low- dimensional and exhibiting certain 
periodicities. Indeed, theoretical investigations of chaotic billiards subject to 
external fields [27], of periodic Lorentz gases [28,29], and of pendulum-like 
differential equations [30,31,32,33,34,35,36] showed that many properties of 
deterministic transport in maps carry over to these more complex chaotic 
dynamical systems. 

In this framework, recently a new feature of deterministic diffusion was discov- 
ered. For simple one-dimensional hyperbolic maps it was shown that the diffu- 
sion coefficient is typically a fractal function of control parameters [37,38,39,40]. 
Subsequently an analogous behavior was detected for other transport coeffi- 
cients [41,42], and in more complicated models [27,28,29]. However, up to now 
the fractality of transport coefficients could be assessed for hyperbolic sys- 
tems only, whereas, to our knowledge, the fractal nature of classical transport 
coefficients in the broad class of nonhyperbolic systems was not discussed. 

On the other hand, studying nonhyperbolic dynamics appears to be more rel- 
evant in order to connect fractal transport coefficients to some known exper- 
iments. Here we think particularly of dissipative systems driven by periodic 
forces such as Josephson junctions in the presence of microwave radiation 
[43,44,45,46,47,48,49], superionic conductors [50,51], and systems exhibiting 
charge- density waves [52] in which certain features of deterministic diffusion 
were already observed experimentally. For these systems the equations of mo- 
tion are typically of the form of some nonlinear pendulum equation. In the 
limiting case of strong dissipation they can be reduced to nonhyperbolic one- 
dimensional time-discrete maps sharing certain symmetries [53,54]. The so- 
called climbing sine map is a well-known example of this class of maps [5,6,9]. 

In this paper we pursue a detailed analysis of the diffusive and dynamical 
properties of the climbing sine map. Particularly, we show that the nonhyper- 
bolicity of this map does not destroy the fractal characteristics of deterministic 
diffusive transport as they were found in hyperbolic systems. On the contrary, 
fractal structures appear for normal diffusive parameters as well as for anoma- 
lous diffusive regions. We argue that higher-order memory effects are crucial to 
understand the origin of these fractal hierarchies in this nonhyperbolic system. 
By using a Green-Kubo formula for diffusion, the dynamical correlations are 
recovered in terms of fractal Takagi-like functions. These functions appear as 
solutions of a generalized integro-differential de Rham-type equation. We fur- 
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thermore show that the distribution of periodic windows exhibiting anomalous 
diffusion forms Devil's staircase like structures as a function of the parameter 
and that the complementary sets of chaotic dynamics have a positive measure 
in parameter space that increases by increasing the parameter value. 

Our paper is organized as follows. In Sec. II we introduce the model. In Sec. 
Ill we explore the coarse functional form of the parameter-dependent diffusion 
coefficient and discuss it in relation to previous results on hyperbolic maps. 
In Sec. IV our analysis is refined revealing complex scenarios of anomalous 
diffusion, which are explained in terms of correlated random walk approxima- 
tions. In Sec. V generalized fractal Takagi functions are constructed for the 
climbing sine map and the connection to the diffusion coefficient is worked 
out. Periodic windows exhibiting anomalous diffusion are studied in detail in 
Sec. VI. We then draw conclusions and discuss our results in the final section. 



2 The climbing sine map 

The one-dimensional climbing sine map is defined as 
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Fig. 1. Illustration of the climbing sine map for the particular parameter value of 
a = 1.189. The dashed line indicates the orbit of a moving particle starting from 
the initial position xq. 

X n+1 = M a (X n ) , M a (X) = X + asin(2vrX) , (1) 

where a G R is a control parameter and X n is the position of a point parti- 
cle at discrete time n. Obviously, M a (X) possesses translation and reflection 
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symmetry, 

M a (X+p) = M a (X)+p, M a (-X) = -M a (X) . (2) 

The periodicity of the map naturally splits the phase space into different boxes, 
(j),p + 1], p e Z, as shown in Fig. f . Eq. (f) as restricted to one box, i.e., on 
a circle, we call the reduced map, 

m a (x) := M a (X) mod f , x := X mod f . (3) 

The probability p n {x)dx to find a particle at a position between x and x + dx 
at time n then evolves according to the continuity equation for the probability 
density p n (x), which is the Frobenius- Perron equation [55] 

Pn+i{x) = J dy p n {y)5{x - m a (y)) . (4) 

The stationary solution of this equation is called the invariant density, which 
we denote by p*(x). 

Due to its nonhyperbolicity, the climbing sine map possesses a rich dynamics 
consisting of chaotic diffusive motion, ballistic dynamics, and localized orbits. 
Under parameter variation these different types of dynamics are highly inter- 
twined resulting in complicated scenarios related to the appearance of periodic 
windows [6,9]. In order to study diffusion we will be interested in parameters 
that are greater than a > = 0.732644... for which the extrema of the map 
exceed the boundaries of each box for the first time indicating the onset of 
diffusive motion. 



3 Coarse structure of the parameter-dependent diffusion coeffi- 
cient 

In this section we explore the relationship between nonlinear maps like the 
climbing sine map and simple piecewise linear maps for which, in contrast to 
the climbing sine map, the diffusion coefficient can be calculated exactly. It 
is well-known that in special cases such different types of maps are linked to 
each other via the concept of conjugacy. Indeed, we show that maps which 
are approximately conjugate to each other exhibit a very similar oscillatory 
behavior in the diffusion coefficient on coarse scales. Our argument refers to 
some existing methods for calculating the diffusion coefficient of piecewise 
linear maps, which we briefly review. We then describe how we numerically 
calculated the complete parameter dependence of the diffusion coefficient for 
the climbing sine map and discuss a first result. 
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3.1 Computing and comparing the diffusion coefficient for approximately con- 
jugate maps 



One speaks of normal deterministic diffusion if the mean square displacement 
of an ensemble of moving particles grows linearly in time. The diffusion coef- 
ficient is then given by the Einstein relation 

D(a)=\im(Xl)/(2n), (5) 

where the brackets denote an ensemble average over the moving particles. 

There exist various efficient numerical as well as, for some system parame- 
ters, analytical methods to exactly compute diffusion coefficients for piecewise 
linear hyperbolic maps, such as transition matrix methods based on Markov 
partitions [37,38,39,40], cycle expansion methods[15, 16, 17,18], and more re- 
cently a very powerful method related to kneading sequences [42]. 

We first restrict our analysis of diffusion in the climbing sine map to parame- 
ters for which there are simple Markov partitions. For one-dimensional maps, a 
partition is a Markov partition if and only if parts of the partition get mapped 
again onto parts of the partition, or onto unions of parts of the partition, see 
Ref. [38] and further references therein. An example of a Markov partition 
consisting of five parts is shown in the inset of Fig. 2. In case of the climbing 
sine map Markov partitions can be constructed simply by forward iteration 
of one of the critical points x c defined by the condition that m a (x c ) = in 
the reduced map. If higher iterations of this point fall onto a periodic orbit 
a Markov partition exits. Indeed, if a Markov partition is known, for piece- 
wise linear maps the diffusion coefficient can often be calculated analytically 
via calculating the invariant measure of the map or via computing the sec- 
ond largest eigenvalue of the Frobenius-Perron operator written in form of a 
transition matrix. 

One can now identify an infinite series of parameter values corresponding to 
a certain type of Markov partition [37,38,39,40]. For parameter values which 
belong to such a Markov partition series the corresponding invariant densities 
p*(x) have a very similar functional form. Note that, in case of nonlinear maps, 
singularities in the invariant density exactly correspond to the iteration of the 
critical point x c [56]. An example of p*(x) for one series of parameter values 
(marked as filled circles) is shown in the inset of Fig. 2. 

By using respective series of Markov partitions piecewise linear maps can be 
related to nonlinear maps. For this purpose let us consider, along with the 
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climbing sine map, (i) the piecewise linear zig-zag map [7,8,16,17] 

!mix n , < x n < 61 

-m 2 (x n - 61) + a, b 1 <x n <b 2 (6) 
mi(x„ - 1) + 1, 6 2 < < 1 

with mi = m 2 = 4a — 1, 61 + b 2 — 1 and ^ = a /mi, and (ii) the nonlinear 
cubic map [6,55] 

3 1 

M a (x n ) = ax 3 n - -ax 2 n + x n (l + -a). (7) 

The definitions of both maps are given on the unit interval. 

In order to compare the diffusion coefficient of these different maps, the pa- 
rameters a were chosen such that the maps all display the same height h, 
defined as the distance between the first iteration of the leftmost critical point 
m a( x c) and the zero bound in the first box (0, 1]. Thus, h — 1 corresponds to 
the onset of diffusion for all three maps. 

For the two simple Markov partition series m a (x c ) = and m a (x c ) = 0.5, 
corresponding to integer and half-integer values of h, respectively, the diffusion 
coefficient of the zig-zag map can be calculated analytically. For integer values 
of h the result reads [8,16] 

_ h(h-l)(2h-l) 
° {h) ~ 3(4h - 1) • (8) 

For half-integer values of h the diffusion coefficient is easily calculated analyt- 
ically, e.g., by transition matrix methods [38,40], to 

In case of the climbing sine map and of the cubic map the diffusion coefficient 
was obtained from computer simulations by evaluating the mean square dis- 
placement Eq. (5) for the same series of Markov partitions. Results are shown 
in Fig. 2. For the climbing sine map some more Markov partition series points 
(alltogether five different series) were included. For all three maps there is a 
very analogous oscillatory behavior of the parameter-dependent diffusion co- 
efficient. These oscillations can be explained in terms of the changes of the 
microscopic dynamics under parameter variation, that is, whenever there is 
a local maximum there is an onset of strong backscattering in the dynamics 
yielding a local decrease of the diffusion coefficient in the parameter, and vice 
versa at local minima [37,38,39,40]. However, the five Markov partition series 
for the climbing sine diffusion coefficient already indicate that there are more 
irregularities on finer scales. For piecewise linear maps, the origin of these 
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irregularities was identified to be the topological instability of the dynamics 
under parameter variation [38,40]. That is, a small deviation of the parameter 
changes the Markov partition and the corresponding invariant density which, 
in turn, is reflected in a change of the value of the diffusion coefficient. Note 
that the dependence of the diffusion coefficient for a single Markov partition 
series appears to be a monotonously increasing function of the parameter [57]. 
Nevertheless, computing D(a) for more and more Markov partitions series will 
reveal more and more irregularities in D(a) thus forming a fractal structure 
[37,38,39,40]. 

Since the climbing sine map shares the same topological features as piecewise 
linear maps in terms of these series of Markov partitions, one may wonder 
whether it is not possible to straightforwardly calculate the diffusion coefficient 
for nonlinear maps from the one of piecewise linear maps by using the concept 
of conjugacy [9,55,58], see also the definition in Appendix A. In fact, it was 
stated by Grossmann and Thomae [9] that the diffusion coefficient is invariant 
under conjugacy, however, without giving a proof. In Appendix A such a 
proof is provided. Unfortunately, conjugacies are explicitly known only in very 
specific cases and for maps acting on the unit interval [55,58]. As soon as 
the map extrema exceed the unit interval, which is reminiscent of the onset 
of diffusive behavior, only some approximate, piecewise conjugacies could be 
constructed in a straightforward way, see Ref. [9] for an example. 

We now apply this reasoning along the lines of conjugacy in order to under- 
stand the similarities between the diffusion coefficient of the three maps as 
displayed in Fig. 2. The functional form of the cubic map can be obtained 
from a Taylor series expansion of sin(x n ) by keeping terms up to third or- 
der thus representing a low-order approximation of the climbing sine map. 
This seems to be reflected in the fact that at any odd integer parameter 
value of h the climbing sine map has an invariant density whose functional 
form is very close to the one of the cubic map at parameter value h — 1, 
p*(x) = 7r _1 (a;(l — x))~ l l 2 . Hence, one may expect that both diffusion coeffi- 
cients are possibly trivially related to each other, however, note the increasing 
deviations between the respective results at larger h. 

For h = 1, the cubic map and the piecewise linear zig-zag map are now in turn 
conjugate to each other [55,58]. However, for ft > 1 we are not aware of the 
existence of any exact conjugacy between zig-zag and cubic map. Still, along 
the lines of Ref. [9] one can at least approximately relate both maps to each 
other via using piecewise conjugacies. This explains why the zig-zag map and 
the climbing sine map display qualitatively the very same oscillatory behavior 
in the diffusion coefficient, somewhat linked by diffusion in the cubic map. 

In summary, by using Markov partitions and by arguing with the concept of 
conjugacy we have shown that the structure of the diffusion coefficient for the 
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nonlinear climbing sine map has much in common with the one of respective 
piecewise linear maps, in the sense of displaying a non-trivial oscillatory pa- 
rameter dependence. However, to use conjugacies in order to exactly calculate 
the diffusion coefficient for nonlinear maps does not appear to be straightfor- 
ward [56], hence in the following we restrict ourselves to alternative methods 
as discussed in the next subsection. 





Fig. 2. Diffusion coefficient at certain Markov partition parameter values for three 
different maps, which are the zig-zag map Eq.(6) (squares), the climbing sine map 
(circles), and the cubic map Eq. (7) (triangles). The values for the zig-zag map 
represent analytical results, see Eqs.(8,9), the remaining values are from computer 
simulations. The lines are guides for the eyes. The inset shows an example of a 
Markov partition for the climbing sine map on the unit interval and the corre- 
sponding invariant density. 



3.2 Complete parameter dependence of the climbing sine diffusion coefficient 



In order to obtain the full parameter dependence for the diffusion coefficient 
of the climbing sine map we numerically evaluated the Green-Kubo formula 
for diffusion in maps [3,4,8,29,38,41] reading 

D n (a) = (j a (x)J:(x))- 1 -(f a (x)). (10) 

Here the angular brackets denote an average over the invariant density of the 
reduced map, (...)= f dxp*(x) . . .. The jump velocity j a (x) is defined by 

j a (x n ) := [X n+1 ] - [X n \ = [M a {x n )\ , (11) 

where the square brackets denote the largest integer less than the argument. 
The sum 



8 




Fig. 3. Diffusion coefficient for the climbing sine map over a large range of pa- 
rameter values. Note the oscillations on large and small scales. The small scale 
fluctuations represent regions of anomalous diffusion where the diffusion coefficient 
either diverges or vanishes. Some of the divergent regions are cut off after a certain 
number of iterations showing plateaus instead of singularities. The data set consists 
of 265005 points. 



k=0 

gives the integer value of the displacement of a particle after n time steps that 
started at some initial position x, and we call it jump velocity function. Eq. 
(10) defines a time-dependent diffusion coefficient which, in case of normal 
diffusion, converges to 

D(a) = lim DJa) . (13) 

n— »oo 

In our simulations for Fig. 3 we truncated J™ (x) after 7 time steps. The invari- 
ant density was obtained by solving the continuity equation for p*(x) Eq. (4) 
with the histogram method of Ref. [1]. Note that both the integer displacement 
and the density are coupled via Eq. (10). Results for D(a) are shown in Fig. 
3 for a large range of parameters demonstrating a highly non-trivial behavior 
of the diffusion coefficient. The large-scale oscillations as predicted from the 
simple Markov partition series (see Fig. 2) are still clearly seen, however, on 
top of this there exist further oscillations on finer scales. These are regions of 
anomalous diffusion that manifest themselves in form of abrupt divergences, 
D(a) — > oo, or by a vanishing diffusion coefficient, D(a) — > 0. 
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4 Simple and correlated random walk approximations 

In this section we study the parameter-dependent diffusion coefficient in more 
detail. Based on the Green-Kubo formula we derive a systematic hierarchy of 
approximations for the diffusion coefficient and show how they can be used to 
understand the complex behavior of this curve in more detail. 

4-1 Asymptotic functional form of the diffusion coefficient on large scales 

We are first interested in understanding the coarse functional form of the 
parameter-dependent diffusion coefficient in the limit of very small and very 
large parameter values. For this purpose we use simple random walk approx- 
imations that are based on the assumption of a complete loss of memory 
between the single jumps. Such an analysis was already performed for hyper- 
bolic piecewise linear maps [38,39]. Here we apply the same reasoning to the 
nonlinear case of the climbing sine map. 

We start in the limit of very small parameter values, i.e., near the onset of 
diffusion. Here we assume that particles make either a step of length one to 
the left or to the right, or just remain in the box. The transition probability 
is then given by integrating over the respective invariant density in the escape 
region. Putting all this information into Eq. (5) yields [6] 

D(a)^p(x c )(2e/a n 2 ) 1 / 2 , (14) 

where e = a — a . Making the additional approximation that p(x) ~ 1 we get 

D(a) ~ 0.525e 1/2 , e < 1. (15) 

The other limiting case concerns values of a > 1. Here the precise value of 
the width of the escape region is much less important than the precise value 
of the step length which is very large, hence by again assuming that p(x) ~ 1 
Eq. (5) can be approximated to [38,39] 

1 r 1 a 2 

D(a)m- dx (M a (x) - xf « — , a > 1. (16) 

2 Jo 4 

These two asymptotic random walk approximations are shown in Fig. 4 as 
bold lines. One can clearly see that there is a dynamical crossover between the 
different functional forms of these two asymptotic regimes. This crossover was 
first observed in piecewise linear maps and appears to be typical for diffusive 
systems exhibiting some spatial periodicity [38,39]. It was lateron also verified 
for diffusion in the periodic Lorentz gas [28] . The coarse functional form of the 
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Fig. 4. Dynamical crossover in the climbing sine map. The two asymptotic regimes 
at small and large parameter values are approximated by simple random walk argu- 
ments yielding the functional forms shown as bold lines. The series of more irregu- 
lar curves corresponds to higher-order correlated random walk approximations with 
n = 0, 1, 2, 3, which give an indication of the exact functional form. The magnifica- 
tion (inset in normal scale) of a small region around a = 3 shows more irregularities 
on a finer scale pointing towards a fractal structure of the diffusion coefficient. 

random walk approximations should be compared to a certain series of higher- 
order approximations based on the Green-Kubo formula Eq. (10), which is 
also shown in the figure. These refined approximations are getting closer and 
closer to the exact functional form, as explained below, and thus give a good 
indication for these exact values. 

4-2 Fine scale of the diffusion coeffcient: anomalous diffusion and bifurca- 
tions 

In the previous subsection we have outlined two simple random walk approx- 
imations for diffusion that do not include any memory effects. However, one 
can do better by systematically evaluating the single terms as contained in the 
series expansion of the Green-Kubo formula Eqs. (10,12). For a simple piece- 
wise linear map and for the periodic Lorentz gas this was done in Ref. [29] and 
provided a simple approach to understand the origin of complex structures in 
the diffusion coefficient on fine scales. 

The basic idea of this approach is as follows: The Green-Kubo formula Eq. 
(10) splits the dynamics into an inter-box dynamics, in terms of integer jumps, 
and into an intra-box dynamics, as represented by the invariant density. We 
first approximate the invariant density in Eq. (10) to p(x) ~ 1 irrespective 
of the fact that it is typically a very complicated function of x and a [6,38]. 
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The resulting approximate diffusion coefficient we label with a superscript 
in Eq. (10), D^(a). The term for n = obviously excludes any higher-order 
correlations and was already worked out in form of the simple random walk 
approximation Eqs. (14)-(16). 

The generalization D^a) , n > 0, which systematically includes more and 
more dynamical correlations, may consequently be denoted as correlated ran- 
dom walk approximation [29]. We now use this expansion to analyze the pa- 
rameter dependence of the diffusion coefficient of the climbing sine map in 
terms of such higher-order correlations. Fig. 4 depicts results for D n (a) at 



Fig. 5. (a) Sequence of correlated random walks D*(a), see Eq. (10) with uniform 
invariant density p{x) ~ 1, for n = 1, • • ■ ,10. Note the quick convergence for normal 
diffusive parameters. The dashed lines define periodic windows, which are the same 
as in Fig.6. The insets (b) and (c) contain two magnifications of (a) in the region 
close to the onset of diffusion for D\ (a) only. They show self-similar behavior on 
smaller and smaller scales. 

n = 0,1,2,3 over a large range of parameters, whereas Fig. 5 (a) presents 
a respective detailed analysis for the region close to the onset of diffusion, 
i.e., for parameters a G (0.732644,1.742726], showing results for D n (a) at 
n — 1, . . . , 10. The series of approximations in Fig. 4 clearly reveals finer and 
finer sequences of oscillations that eventually converge to a fractal structure, 
as is particularly shown in the inset of this figure. However, here the order of 
the expansion is not large enough to identify parameter regions of anomalous 
diffusion. These regions can be better seen in Fig. 5, where three different 
cases of parameter regions can be distinguished: (i) regions with quick conver- 
gence of this approximation corresponding to normal diffusion (ii) divergence 
of -D^(a) corresponding to ballistic motion, in agreement with D(a) — > oo, and 
(iii) localized dynamics where D^(a) alternates in n between two solutions, 
with D(a) — > for the exact diffusion coefficient. This oscillation points to 
the dynamical origin of localization in terms of certain period-two orbits. That 
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these approximate solutions are non-zero is due to the fact that the invariant 
density was set equal to one. The dashed lines in Fig. 5 indicate the largest 
regions of anomalous diffusion. The approximate diffusion coefficient D{ Q (a) 
of this figure is compared to the "numerically exact" one in Fig. 6. Here "nu- 
merically exact" we wish to be understood in the sense that no further ad 
hoc-approximations are involved, i.e., we evaluated the Green-Kubo formula 
according to the numerical method described in Sec. III. by truncating it after 
20 time steps. This comparison shows that in case of normal diffusion our 
approximation nicely reproduces the irregularities in the non-approximated 
diffusion coefficient. Like the inset of Fig. 4, the magnifications in Fig.5 give 
clear evidence for a self-similar structure of the diffusion coefficient. These 
results thus show that dealing with correlated jumps only yields a qualitative 
and to quite some extent even quantitative understanding of the regions of 
normal and anomalous diffusion in the climbing sine map. 

The impact of specific features of the microscopic dynamics on the diffusion 
coefficient is nicely elucidated by comparing the bifurcation diagram of the 
reduced climbing sine map Eq. (3) with the numerically exact diffusion coef- 
ficient, see Fig. 6. As one can see in the upper panel of Fig. 6, the bifurcation 
diagram consists of (infinitely) many periodic windows. Whenever there is a 
window the dynamics of Eq. (1) is either ballistic or localized [6,8]. Fig. 6 
demonstrates the strong impact of this bifurcation scenario on the diffusion 
coefficient. For localized dynamics, orbits are confined within some finite in- 
terval in phase space implying subdiffusive behavior for which the diffusion 
coefficient vanishes, whereas for ballistic motion particles propagate superdif- 
fusively with a diverging diffusion coefficient being proportional to n. Only 
for normal diffusion D(a) is nonzero, finite and the limit in Eq. (13) exists. 
At the boundaries of each periodic window the diffusion coefficient is related 
to intermittent-like transient behavior eventually resulting in normal diffusion 
with D(a) ~ a^) [5,6,9,10]. 



5 Diffusion coefficient in terms of fractal generalized Takagi func- 
tions 

In this section we further analyze the dynamical origin of the different struc- 
tures in the parameter-dependent diffusion coefficient by constructing objects 
called fractal generalized Takagi functions. These functions somewhat resem- 
ble usual Takagi functions but, as will be shown, they fulfill a more complicated 
type of functional recursion relations than standard de Rham-type equations. 
Interestingly, Takagi functions were known to mathematicians since about a 
hundred years [59,60,61], however, in the field of chaotic transport they were 
appreciated by physicists only very recently [3,4,38,62,63]. 
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Fig. 6. Upper panel: bifurcation diagram for the climbing sine map. Lower panel: 
diffusion coefficient computed from simulations as a function of the control parame- 
ter a in comparison with the correlated random walk approximation D\ (a) (dots). 
The dashed vertical lines connect regions of anomalous diffusion, D(a) — > oo or 
D(a) — > 0, with 'ballistic' and 'localized' windows, respectively, of the bifurcation 
diagram. 

We first show how to construct fractal generalized Takagi functions and study 
their properties with respect to the three different types of dynamics in the 
climbing sine map. We then relate these objects directly to the diffusion coef- 
ficient . 

From the definition of the time-dependent jump velocity function there follows 
the recursion relation [37,38] 

J2(x)=j a (x) + J2- 1 (m a (x)). (17) 

Since the time-dependent jump velocity function J™ is getting extremely com- 
plicated after some time steps, we introduce the more well-behaved function 

:= f J:(z) dz, T a «(0) = T a «(l) = 0. (18) 
J o 

Integration of Eq. (17) yields the recursive functional equation 

T:(x) = t a (x) + ~~tt~\ Tr\m a (x)) - I(x) (19) 

rn a\ X ) 

with the integral term 

I(x) = / dzg" a {z)T:-\z) , (20) 
Jo 
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Fig. 7. (a): Generalized fractal Takagi functions for the diffusive climbing sine map 
with parameters a = 1.2397 (upper curve) and a = 1.7427 (lower curve), (b), (c): An 
example of nonconverging iterations of the generalized fractal Takagi functions for 
the climbing sine map with parameters corresponding to (b) ballistic dynamics at 
a = 1.0 and to (c) localized dynamics at a = 1.5, both for the time steps n = 5, 6, 7. 
Note the divergence of the iterations in (b) and the alternation between two states 
in (c). 

function t a (x) is given by 

t a (x) = / dzj a (z) = Xj a (x) + c(x) , (21) 
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where c(x) is defined to be constant on each subinterval where the jump ve- 
locity j a {x) has a given value. This constant is fixed by the condition for t a (x) 
to be continuous on the unit interval supplemented by 

ta(0) = t a (l) EE 0. (22) 

The generalized Takagi function is now defined in the long-time limit of Eq. 

(19), 

T a (x) = Vml?(x) . (23) 

For piecewise linear hyperbolic maps I(x) in Eq. (20) simply disappears, and 
the derivative in front of the second term of Eq. (19) reduces to the local slope 
of the map thus recovering ordinary de Rham-type equations [3,4,38,41,42]. 
It should be noted that for smooth nonlinear maps like the climbing sine 
map the reduced map m a (x) is generally not invertible. In order to define a 
local inverse of m a (x), we split the unit interval into sub intervals on which this 
function is piecewise invertible. Thus Eq. (19) should be understood as a series 
of equations where each part is defined for a respective piecewise invertible 
part of m a (x). The detailed derivation of Eq. (19) is given in Appendix B. 

It is not known to us how to directly solve this generalized de Rham-type 
integro-differential equation for the climbing sine map, however, solutions can 
alternatively be constructed from Eqs. (18), (21), (22) starting from simula- 
tions. Results are shown in Fig. 7. For normal diffusive parameters the limit 
of Eq. (23) exists and the respective curves are fractal on the whole unit inter- 
val somewhat resembling standard Takagi functions [3,4,38,62,63]. However, 
in case of periodic windows T£(x) either diverges due to ballistic flights, or 
it oscillates indicating localization. Interestingly, in these functions the corre- 
sponding attracting sets appear in form of smooth, non-fractal regions on fine 
scales as marked by the dashed lines in Fig. 7. 

The diffusion coefficient can now be formulated in terms of these fractal func- 
tions by carrying out the integrations contained in Eq. (10). For simplicity we 
restrict ourselves to the parameter region of a G (0.732644, 1.742726] in which 
the respective solution reads 

D(a) = 2 [T a (x 2 )p(x 2 ) - T a (x 1 )p(x l )] - D p (a), (24) 

where Xi, % — 1,2, is defined by [M (xj)] := 1, and Dq(o) := dxp(x). Our 
previous approximation D^(a) with p{x) ~ 1 is recovered from this equation 
in form of 

Dl{a) = 2 [T a (x 2 ) - T a ( Xl )} - D^a) . (25) 

Hence, Eqs. (24), (25) explicitly relate the generalized fractal Takagi functions 
shown in Fig. 7 to the fractals of Fig. 6. 
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6 Periodic windows 



One of the most important problems regarding periodic windows in maps 
remains the question of their total measure. Much understanding has been 
achieved for one-dimensional unimodal maps [64,65,66,67,68,69,70]. Based on 
the Sharkovskii theorem about the ordering of periodic orbits [64], Metropolis, 
Stein and Stein organized periodic windows in universal symbolic sequences 
(U-sequences) such that the sequence of next order is uniquely determined by 
the previous one [65]. Later Jacobson came up with the proof that chaotic 
parameter values in one-dimensional unimodal maps with a single maximum 
do have positive measure [66]. Related numerical studies were made by Farmer 
[69]. Furthermore, it was shown that periodic windows in such a map form 
so-called fat fractal Cantor-like sets with positive measure. 

However, for diffusive maps on a line, apart from the preliminary studies of 
Refs.[6,8], nothing appears to be known. On the other hand, as was exemplified 
in Sec. IV there is an intimate relation between the irregular structures of the 
diffusion coefficient and the occurrence of periodic windows. Hence, in this 
section we investigate the periodic windows for the climbing sine map in full 
detail. 

Due to the spatial extension of our model a new type of periodic motion, 
which is not present in unimodal maps, exists, which is that particles move 
on average in one direction. These ballistic solutions are born through tangent 
bifurcations, further undergo a Feigenbaum-type scenario and die at crises 
points [6,8]. Localized solutions occurr at even periods only and start with 
tangent bifurcations followed by a symmetry breaking at slope-type bifurca- 
tion points [6,8]. In this case the bifurcation scenario is much more complex. 
Obviously, periodic windows are related to certain periodic orbits, thus there 
are infinitely many of them, and they are believed to be dense in the parameter 
space [2]. 

By dividing the parameter line into subsets labeled by the integer value of the 
map maximum on the unit interval, [M a (X cr )] = j, j G Z, we computed all 
windows of a certain period p in a certain subset. The numerical procedure 
which was used for these computations is outlined in Appendix C. Let Tj 3 
denote the total measure of all period p-windows in a subset j and let Sj be 
the partial sum of Tf defined by = £f=i Tj. In Fig. 8 log 7? is plotted as 
a function of the period for the three first subsets j = 0, 1,2. The measures 
corresponding to windows with localized orbits are shown in Fig. 8 as pluses. 
Is it clear that they make the major contribution to the total measure for even 
periods hence explaining the origin of the pronounced oscillatory behavior of 

T3 
p' 
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Fig. 8. Upper panel: The total measure Tj of all period p-windows (lines with 
symbols) for the first three parameter intervals (from top to bottom) as defined in 
the text. The dotted lines represent exponential fits; for the parameters see Table 1. 
The measures corresponding to windows with localized orbits are shown as pluses. 
Lower panel: The partial sum 5j for all periodic windows at a certain period p. The 
dashed curves represent approximations as calculated from Eq. (27), the straight 
lines are their limiting values at p — > oo. 

Table 1 

Fit parameters for the exponential decrease of the measure at even and odd periods 
for the first three subsets of the map control parameter labeled by j. 



j 


„even 
a j 


^even 


n odd 
U 3 


yadd 





0.284 


0.61 ±0.01 


0.012 


0.62 


1 


0.224 


0.87 ±0.01 


0.007 


0.89 


2 


0.242 


1.01 ±0.02 


0.006 


0.98 



For odd periods all windows are due to ballistic orbits. Thus, we more carefully 
distinguish between even and odd periods. The dotted lines in Fig. 8 represent 
exponential fits to the functional dependence of the measure at even and odd 
periods according to a e ^ en exp(— bj ven p) and a° dd exp(— b° dd p), where j stands 
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Fig. 9. Upper panel: Total measure Tj for period one-windows as a function of the 
control parameter interval labeled by j. The slope of the line appears to be —1 (see 
text). Lower panel: the same as on the upper panel but here for different window 
periodicities p, and for smaller j. The labels at the single graphs give the value of 
p. 

for the box number j = 0, 1, 2, as denned above. The fit parameters a,j, bj are 
given in Table 1. From Fig. 8 and Table 1 one can conclude that the total 
measures at even and odd periods at a certain label j decrease approximately 
with the same rate, bj Ven ~ b° dd ~ bj. The exponential decrease of Tj suggests 
that the measure of chaotic solutions in each box, which is complementary to 
the measure of periodic windows, is indeed positive. Based on the information 
of Tj, it is straightforward to approximate the total measure of all the periodic 
windows in the jth box by 

oo ' oo " „even _i_ „odd p —bj 

S i = T}+Y,Tj + -£^^T}+ ' (26) 

i=2 i=3 

and to approximate the parameter dependence of the partial sum to 

~ T) e- b > ( p - 5 ) + Sj ( 1 - e~ b i {p - s } ) . (27) 
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Table 2 

The Lebesque measure Aj, the total measure of periodic windows Sj, and the com- 
plementary measure of chaotic solutions Cj for the first three subsets of the map 
control parameter labeled by j. 



3 


A? 




Cj 





1.01008 


0.226 ± 0.002 


0.783 ± 0.002 


1 


1.00265 


0.103 ±0.002 


0.898 ± 0.002 


2 


1.00123 


0.068 ± 0.002 


0.932 ±0.002 



Sums with one or two primes go only over even or odd terms, respectively. 
In the lower panel of Fig. 8 results for Sj are shown as horizontal lines, the 
dashed lines in the upper panel (without symbols) are the approximations for 
Sj according to Eq. (27). 

The values for the measure of all periodic windows in the jth box, Sj, and the 
measure of the chaotic solutions Cj = Aj — Sj, where Aj is the total measure 
of the j box, are listed in Table 2 for the first three subsets of the control 
parameter. 




J T I I I I I I I 

0.8 1 1.2 1.4 1.6 

a 



Fig. 10. Devil's staircase like structure formed by the distribution of periodic win- 
dows as a function of the control parameter. N is the integrated number of period 
six-windows. The inset shows a blowup of the initial region. 

The main message of Table 2 is that the measure of periodic windows is 
different for a different subset j and obviously decreases by increasing the 
control parameter. Correspondingly, the measure of chaotic solutions increases 
as the parameter of the system is getting larger. In order to make this more 
quantitative we consider the dependence of the measure of period p-windows 
as a function of the box index j. First we check period one-windows. Since in 
each box there exists only one window of this period, we are able to go to up 
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to j = 20. In the upper panel of Fig. 9 the logarithm of Tj is plotted against 
log j. We find that the slope of this function almost exactly equal to —1. The 
behavior of Tj for different p is shown in the lower panel of Fig. 9. The slope 
of the line for period two is also — 1, for period three and four it is —2, and 
for period five it is —3. Fig. 9 shows that even and odd periods decrease with 
respectively different laws, where the decay rate appears to be precisely given 
by the periodicity of the windows according to l/j p ^ 2 for even periods and 
1/j'(p +1 )/ 2 for odd periods. 

In order to analyze the structure of the regions of anomalous diffusion in the 
parameter space, we sum up the number of period six-windows as a function of 
the parameter, that is, the total number is increased by one for any parameter 
value at which a new period six-window appears. This sum forms a Devil's 
staircase like structure in parameter space indicating an underlying Cantor 
set like distribution for the corresponding anomalous diffusive region, see Fig. 
9. Since the Lebesque measure of periodic windows is positive, this set must 
be a fat fractal [68]. Its self-similar structure can quantitatively be assessed 
by computing the so-called fatness exponent. Following [69], let h(e) be the 
total measure of all periodic windows whose width is greater than or equal 
to e. Define the coarse-grained measure as fx(e) = A — h(e), where A is 
the total measure of a box related to the control parameter. For quadratic 
maps on the interval, it was conjectured and confirmed numerically that fi(e) 
asymptotically scales as a power law in the limit of e — > 0, 

fi(e) « /i(0) + Ae p , (28) 

where /x(0) is the measure of chaotic parameters. f3 was called the fatness expo- 
nent. For quadratic maps it was found to be f3 ~ 0.45. Since our map belongs to 
the same universality class as considered in Farmer's case, namely the map has 
a single quadratic maximum, one may expect that (3 will have the same value. 
For the climbing sine map a double-logarithmic plot of 5fi(e) = fi(e) — /x(0) is 
shown in Fig. 11 for the first three boxes, j = 0, 1, 2. In all cases the bold lines 
have the slope 0.45 with errors of 0.03, 0.04, 0.05, respectively, which seems to 
be in agreement with Farmer's conjecture about the universality of (3. However, 
apart from this coarse linear behavior one can see an interesting oscillatory 
behavior in 5/i with respect to e. This fine structure can be explained with 
respect to the histogram distribution functions / of the window sizes d, which 
are plotted in Fig. 11 in form of dashed lines. Somewhat surprisingly, the peri- 
odic windows are not distributed uniformly or smoothly with respect to their 
size but form certain peaks, in which preferably windows of certain periods are 
grouped together. This non-uniformity is clearly reflected in the oscillations of 
5/i(e). Moreover, we find that the size distributions of periodic windows have 
a fine structure that appears to resemble a fractal function. Some evidence 
for this property is given in the inset of Fig. 11 (c), which shows the size of 
every window of period five in the subset j = 2 as a function of its appearance 
with respect to the map control parameter a, i.e., not the parameter itself 
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Fig. 11. Difference Sfi = fi(e) — /u(0), where is the total measure of windows that 
are smaller than e. Results are plotted for the first three subsets j = 0, 1, 2 of the 
map control parameter, from top to bottom. In all cases, the slope of the solid lines 
is approximately 0.45. The dashed lines in each graph represent the corresponding 
coarse-grained window distribution functions. In the inset of Fig. (c) the measure 
of any single period five-window in the subset j = 2 is shown with respect to an 
integer label that accounts for the ordering according to the map control parameter. 



is plotted but just an integer running index iV is given instead. Particularly 
the height of the peaks is important, and one can clearly see a complicated 
hierarchy of different peaks which are reminiscent of the fine structure in the 
corresponding distribution function shown in Fig. 11 (c). 
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7 Conclusions 



In this paper we have performed a detailed analysis of the parameter depen- 
dence of the diffusion coefficient in a nonhyperbolic dynamical system. The 
climbing sine map has been chosen as a paradigmatic example of such a sys- 
tem. We have shown that, on a coarse scale, there are certain analogies between 
the parameter-dependent diffusion coefficient of this map and the ones in sim- 
ple hyperbolic piecewise linear maps, such as the existence of an oscillatory 
structure, and the existence of asymptotic functional forms as derived from 
simple random walk models. However, in contrast to hyperbolic maps showing 
normal diffusion only, in the nonhyperbolic climbing sine map fractal struc- 
tures appear for both normal and anomalous diffusive regions of the diffusion 
coefficient. An understanding of the origin of these fractal structures was given 
in terms of dynamical correlations starting from the Green-Kubo formula for 
diffusion. We furthermore related these irregularities in the diffusion coeffi- 
cient more microscopically to different characteristics in corresponding fractal 
generalized Takagi functions. For this purpose we derived a new functional 
recursion relation that defines these fractal forms and generalizes ordinary 
de Rham-type equations. Our analysis was completed by extensive numeri- 
cal studies of the periodic windows of the climbing sine map showing that 
both the periodic and the chaotic parameter regions have positive measures 
in the parameter space. However, these measures are themselves parameter- 
dependent, and by increasing the parameter we found that the chaotic regions 
occupy larger and larger measures. We finally provided evidence that these 
different sets form fat fractals on the parameter axis. 

In conclusion, we wish to remark that the climbing sine map is of the same 
functional form as the respective nonlinear equation in the two-dimensional 
standard map, which is considered to be a standard model for many physical 
Hamiltonian dynamical systems. Indeed, both models are motivated by the 
driven nonlinear pendulum, both are strongly nonhyperbolic, and though the 
standard map is area-preserving it too exhibits a highly irregular parameter- 
dependent diffusion coefficient. Understanding the origin of these irregularities 
was the subject of intensive research [2,19], however, so far the complexity of 
this system did not enable to reveal its possibly fractal nature. A suitably 
adapted version of our approach to nonhyperbolic diffusive dynamics as pre- 
sented in this paper may enable to make some progress in this direction. 

Another interesting problem is to possibly further exploit the concept of con- 
jugacy between nonlinear and piecewise linear maps, as explained in Sec. Ill, 
in order to exactly calculate diffusion coefficients for nonlinear maps. A very 
promising approach in this direction was presented in Ref . [56] . Based on these 
techniques we are planning to perform a spectral analysis of the Frobenius- 
Perron operator governing the probability density of the diffusive climbing sine 
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map. Combining such an analysis with the Takagi function approach outlined 
here may lead to a general theory of nonhyperbolic transport. 

It would furthermore be important to check out the applicability of periodic 
orbit theory for computing the parameter-dependent diffusion coefficient of 
the climbing sine map, which may provide an alternative method [15]. An- 
other promising direction of future research concerns establishing crosslinks 
between our work and the realm of strange kinetics and stochastic modeling 
as described in Refs. [24,25,26], e.g., by trying to apply continuous time ran- 
dom walk techniques to more complicated chaotic models exhibiting fractal 
diffusion coefficients such as the climbing sine map. 

We finally emphasize the importance to look for possibly fractal transport 
coefficients in experiments. A very promising candidate appears to be the 
phase dynamics in SQUID's, which was very recently analyzed theoretically 
[48] and studied experimentally [49]. 



A Diffusion coefficients of two conjugate maps 

In this Appendix we give a proof of the statement of Grossmann and Thomae 
[9] that two diffusive maps which are conjugate to each other have the same 
diffusion coefficient. 

Two diffusive maps F : / — > / and G : J — > J are called conjugate [58,9,55] if 
there exists a map H : / — > J such that F(x) = H(G(H~ 1 (x))). Let us assume 
in the following that the conjugation function H is sufficiently smooth. Let 
the invariant densities of the corresponding reduced (mod 1) maps be p(x) for 
F(x) mod 1 and p(y) for G(y) mod 1; then it is, according to conservation 
of probability, p(x) = \(H^ 1 (x))'\ p(H~ 1 (x)). The diffusion coefficients of the 
maps F(x) and G(y) we denote by D F and D G , respectively. Without loss of 
generality let us furthermore assume that the maxima of both maps are in the 
interval [1,2]. 

We now start with the Green-Kubo formula written in correlated random walk 
terms as 

D F = -J [F(x)] 2 ~p(x)dx + J [F(x)} ■ B(x) p(x)dx, (A.l) 



where 

B(x) = [F{F(x)}] + ••• + [F ({F ({... ({F(x)}) ...})})] + ... , (A.2) 
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or shortly 



where 



D F = d Fo + d Fl + d F2 + ... (A. 3) 



i r 



o 

i 



d Fl = \j [F{x)\ [F{F(x)}\ p(x) dx, (A.5) 
o 

and so on. Focusing on the first term, one can rewrite this expression using 
the symmetry of the map to 

1 x2 

dp ° = 2 / [F{X)]2 P{X) dX = J ~ P ^ X) dX ' (A ' 6) 

xl 

where x±,x 2 defines an escape region. For the conjugate map G(y) the respec- 
tive term reads 

1 2/2 

**> = 2 / My)? M d y = J M d y> ( A - 7 ) 

yl 

where 2/1,2/2 is the corresponding escape region for G. Note that the escape 
regions (x\,X2) and (2/1,2/2) are not the same, however, it is straightforward to 
show that H(xi) — j/j , % — 1,2, that is, the topology of both maps is conserved 
such that the two escape regions are mapped onto each other under conjugacy. 

Taking into account the conservation of probability mentioned before one im- 
mediately gets 

d Fo = d Go . (A. 8) 

All other terms d Fl , d F2 , ... and dc 11 dc 2 , ■■■ have the form 

d {F>G)i =A Jv(z)dz, i = 1,2, ... (A.9) 

& esc 

where A is a constant, 5 esc is the respective escape region and u(z) dz is the 
corresponding invariant measure. Thus, the same argument can be applied to 
show that d Fi = do v (i = 1, 2, ...). Combining all results we arrive at 

D F = Dq . (A. 10) 
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B Recursion relation for generalized Takagi functions 



We start with the recursion relation for the jump velocity function Eq. (17), 

J:(x)=j a (x) + J:~ 1 (m a (x)), (B.l) 
by recalling the definition of the generalized Takagi function Eq. (18), 

T:(x) : = f J:(z) dz, 7?(0) = T:(l) = 0, (B.2) 
Jo 

or differently 

= ^(x). (B.3) 




Fig. B.l. Illustration of the construction of the inverse function of the climbing sine 
map for the parameter value a = 1.189. Piecewise invertible branches are labeled 
by the integer numbers i = 1, ■ ■ ■ , 7. 

We have to integrate Eq. (B.l), 

jdy J:(y) = jdyj a (y) + j dy J n a -\m a (x)). (B.4) 



By using of Eq. (B.2) we get 

T:(x)=t a (x) + I(x), xe(0,l], (B.5) 

where 

X 

I(x) = Jdy J n a -\m a {y)), x G (0, 1]. (B.6) 
o 

Without loss of generality let us assume that the maximum of the map is in 
the interval [1,2]. 
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Depending on x the integral in Eq. (B.6) can be decomposed into 

X 

h(x) = Jdy J2-\m a {y)), x G (0,^]; (B.7) 
o 

XI X 

h{x) = Jdy J2-\m a {y)) + J dy J2~\m a (y)), x G (x 1 ,x 2 }; 

xi 

I 3 (x) = x G {x 2 ,x 3 ]; ... I 6 (x) = x G (x 5 ,x e ]; 

X\ X 

h{x) = J dy J2-\m a {y)) + ■■■ + J dy J^ l (m a (y)), x G (x 6 ,x 7 ]. 

x 6 

Each integral in Eq. (B.7) now contains only one piecewise invertible branch 
of the reduced map m l a (x) as shown in Fig. B.l. Here, the piecewise invertible 
branches of the reduced map are labeled by integers, and the corresponding 
branches of the inverse function g l a (y) have the same indices i — 1, • • • , 7. Since 
all integrals in Eq. (B.7) have the same form (only the inverse parts of the 
reduced map are different), we restrict ourselves to the integral 

X 

I(x) =Jdy Jr\K{y)). (B.8) 
o 

Making the change of variables z = m l a (y) and using the definition of the 
generalized Takagi function Eq. (B.3) we get 

m*,(x) 

J dz(gi(z))'-Tr\z). (B.9) 
o 

Using integration by parts we arrive at 

m* (x) 

I(x) = (gi(z))' ■ T:~\z)\Z l{x) - / dz (gl(z))" • T" -1 (z). (B.10) 

o 

Now recall that according to Eq. (B.2) it is T™ -1 (m l a (x j)) = 0, where the Xj 
define the boundaries of the piecewise invertible parts of m a (x), see Fig. B.l, 
and that (g l a (z)) | z=m »( x ) = l/(m l a (x)) . Thus, by formally defining the inverse 
function g a (x) as consisting of all branches i — 1, . . . , 7, we can finally write 
Eq. (B.5) in the form 

T:(x) = t a (x) + -1- T:-\m a (x)) - I(x) (B.ll) 

' a a\ x ) 
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with the integral term 

I(x) = / dzg" a {z)T:~\z). (B.12) 
J o 

C Numerical procedure for calculating the measure of the periodic 
windows 

The parameter values at an which correspond to the tangent bifurcations of 
the p-periodic windows were found by solving the two coupled transcendental 
equations 

dmWjx)/dx = 1, mg» - x = 0, (C.l) 

where m^(x) denotes the p-times iterated reduced map. This corresponds to 
the situation where mffl(x) touches the bisector. Somewhat after a tangent 
bifurcation one will unavoidably find a situation where a critical point x c , 
which corresponds to an extremum of (x) , crosses the bisector. When this 
critical point is exactly located on the diagonal, the reduced map or its higher 
iterations have a fixed point and there exists a specific Markov partition on the 
interval [37,38]. The periodic orbit generated by the corresponding parameter 
value a ss is superstable, 

™£?> c ) - z c = 0. (C.2) 

By further increasing the parameter value up to a cr a crisis takes place, and 
this again corresponds to the existence of a certain Markov partition. 

Based on this scenario, the full numerical procedure which was used for cal- 
culating the measure of periodic windows is as follows: The values of a ss cor- 
responding to superstable solutions were first calculated by a combination of 
bisection with the Newton method. The parameters for the tangent bifurca- 
tions could then usually be found by the modified two-dimensional Newton 
method [69]. However, the highly discontinuous nature of m^(x) made its 
implementation very inefficient. Instead, starting in the vicinity of each a ss 
we again combined the one-dimensional Newton and bisection methods. This 
ensured that no windows were missed. Finally, the parameter values corre- 
sponding to crisis points a cr , which are also defined by Markov partitions, can 
be found by solving respective equations that are formally analogous to Eq. 
(C.2). 
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